* Macros for calculating job accessibilty by residence Census tract;
* Produced by Mark Kutzbach;

* all encompassing macro to be run in primary program;
%macro tract_access(otmyear);

* so far, geo is tract;
/* run geo frame and dist ----------------------------------------------------------------------- */
%if (&run_geoframe_dist.=1) %then %do;
/* create frame of geo pairs to use for city and calculate distances */
    * identify geo pairs;
    data trt_o (keep=stfid trctlat trctlon) 
         trt_d (keep=stfid trctlat trctlon);
        set GEO.intpts_trct (rename=(TRCT=stfid)
                            where=(input(substr(stfid,1,2),2.) = &stcode. 
                            and substr(stfid,3,3) in (&cty.)
                            and substr(stfid,6,6) ne '000000')
                            );
    run;
    * create a match of all possible combinations of geo;
    proc sql;
        create table trt_pairs as
            select orig.stfid as o_stfid, orig.trctlat as o_trt_lat, orig.trctlon as o_trt_lon,
                   dest.stfid as d_stfid, dest.trctlat as d_trt_lat, dest.trctlon as d_trt_lon
                from
                    trt_o as orig,
                    trt_d as dest;
    quit;
    * calculate distances;
    %let pi=3.14159265359;
    data trt_pairs (drop=o_trt_lat d_trt_lat o_trt_lon d_trt_lon);
        set trt_pairs;
        dist=3963.1676 * 2 * arsin(sqrt(
                             (sin((o_trt_lat-d_trt_lat)*&pi./360))**2 +
                             cos(o_trt_lat*&pi./180)*cos(d_trt_lat*&pi./180) *
                             (sin((o_trt_lon-d_trt_lon)*&pi./360))**2
                             ));
        if (dist ne .);
    run;
    * sort for next steps;
    proc sort data=trt_pairs out=OUTDATA.trt_pairs_&met.;
        by d_stfid o_stfid;
    run;
%end; * end frame and distance;

/* run proximity measures ---------------------------------------------------------------------- */
%if (&run_proxmeas.=1) %then %do;
    * auto;
    %if ("&use_prox_time_mpo_auto."="1") %then %do;
    proc sort data=PROX.time_mpo_auto_&met. (rename=(time_mpo_auto=auto
                                      work_stfid=d_stfid 
                                      home_stfid=o_stfid))
              out=time_mpo_auto_aug;
        by d_stfid o_stfid;
    run;
    data time_mpo_auto_aug;
        set time_mpo_auto_aug;
        by d_stfid o_stfid;
        if (first.o_stfid);
        count=1;
    run;
    proc summary data=time_mpo_auto_aug nway;
         class d_stfid;
         var count;
         output out=pairs_obs (drop=_TYPE_ _FREQ_) n(count)=pair_obs_auto;
    run;
    data time_mpo_auto_aug;
        merge time_mpo_auto_aug (drop=count)
              pairs_obs;
        by d_stfid;
    run;
    proc sort data=time_mpo_auto_aug;
        by d_stfid o_stfid;
    run;
    %end;
    * transit;
    %if ("&use_prox_time_mpo_tran."="1") %then %do;
    proc sort data=PROX.time_mpo_tran_&met. (rename=(time_mpo_tran=tran
                                      work_stfid=d_stfid 
                                      home_stfid=o_stfid))
              out=time_mpo_tran_aug;
        by d_stfid o_stfid;
    run;
    data time_mpo_tran_aug;
        set time_mpo_tran_aug;
        by d_stfid o_stfid;
        if (first.o_stfid);
        count=1;
    run;
    proc summary data=time_mpo_tran_aug nway;
         class d_stfid;
         var count;
         output out=pairs_obs (drop=_TYPE_ _FREQ_) n(count)=pair_obs_auto;
    run;
    data time_mpo_tran_aug;
        merge time_mpo_tran_aug (drop=count)
              pairs_obs;
        by d_stfid;
    run;
    proc sort data=time_mpo_tran_aug;
        by d_stfid o_stfid;
    run;
    %end;                        
                                            
    * merge in and create desired supplemental proximity measures by pair (not used);
    %macro loadprox(proxtype);
        %let beginpart=use_prox_;
        %if ("&&beginpart.&proxtype.."="1") %then %do;
            PROX.&proxtype._&met. (in=in_&proxtype. 
                                   rename=(work_stfid=d_stfid home_stfid=o_stfid))
        %end;
    %mend loadprox;
    * merge in data should be sorted by work geo home geo;
    data OUTDATA.proxmeas_&met.;
        merge OUTDATA.trt_pairs_&met. (in=in_trt)
            %if ("&use_prox_time_ctpp."="1") %then %do;
            PROX.time_mpo_auto_&met. (rename=(work_stfid=d_stfid home_stfid=o_stfid))
            %end;
            %if ("&use_prox_time_mpo_auto."="1") %then %do;
            time_mpo_auto_aug (in=in_auto)
            %end;
            %if ("&use_prox_time_mpo_tran."="1") %then %do;
            time_mpo_tran_aug (in=in_tran)
            %end;
            ;
        by d_stfid o_stfid;
        * only retain those in geo frame;
        if in_trt;
        * adjustments to auto travel time;
        %if ("&use_prox_time_mpo_auto."="1") %then %do;
        * within tract;
        if (d_stfid=o_stfid) then auto=0;
        * no data available;
        predict_auto=0;
        if ((d_stfid ne o_stfid) and ((auto=.) or (auto le 0))) then do;
            predict_auto=1;
            auto=dist*3;
        end;
        * could assume some time of accessing car (previously added 5 minutes);
        auto=auto;
        * do not ever make worse than walking;
        if (auto>(dist*20)) then do;
            auto=dist*20;
            predict_auto=2;
        end;
        * speed;
        if (d_stfid ne o_stfid) then speed_auto=(dist/auto)*60;
        %end;
        * adjustments to transit travel time;
        %if ("&use_prox_time_mpo_tran."="1") %then %do;
        * within tract;
        if (d_stfid=o_stfid) then tran=0;
        * no data available;
        predict_tran=0;
        if ((d_stfid ne o_stfid) and ((tran=.) or (tran le 0))) then do;
            predict_tran=1;
            tran=dist*20;
        end;
        * assume some minimum time (often already done);
        if (d_stfid ne o_stfid) then tran=max(10,tran)-5;
        * do not ever make worse than walking;
        if (tran>(dist*20)) then do;
            tran=dist*20;
            predict_tran=2;
        end;
        * speed;
        if (d_stfid ne o_stfid) then speed_tran=(dist/tran)*60;
        * best travel time between auto and transit;
        %if ("&use_prox_time_mpo_auto."="1") %then %do;
        better_tran=0;
        best=auto;
        speed_best=speed_auto;
        if (tran le auto) then do;
            better_tran=1;
            best=tran;
            speed_best=speed_tran;
        end;
        %end;
        %end;   
    run;

    * evaluation of mpo proximity;
    %if ("&evaluate_prox."="1") %then %do;
    
    proc means data=OUTDATA.proxmeas_&met. (where=(d_stfid ne o_stfid));
        title "summary of proximity for &met.";
        var dist 
        %if ("&use_prox_time_mpo_auto."="1") %then %do;
        predict_auto auto speed_auto 
        %end;
        %if ("&use_prox_time_mpo_tran."="1") %then %do;
        predict_tran tran speed_tran
        %if ("&use_prox_time_mpo_auto."="1") %then %do;
        better_tran best speed_best;
        %end;
        %end;
        ;
    run;
        
    %if ("&use_prox_time_mpo_auto."="1") %then %do;
    proc freq data=OUTDATA.proxmeas_&met. (where=(d_stfid ne o_stfid));
        title "imputation rate of auto for &met.";
        tables predict_auto;
    run;
    proc means data=OUTDATA.proxmeas_&met. 
        (where=((predict_auto=0) and (d_stfid ne o_stfid)));
        title "summary of reported auto travel time for &met.";
        var dist auto speed_auto
        %if ("&use_prox_time_mpo_tran."="1") %then %do;
        tran speed_tran
        better_tran best speed_best;
        %end;
        ;
    run;
    %end;
        
    %if ("&use_prox_time_mpo_tran."="1") %then %do;
    proc freq data=OUTDATA.proxmeas_&met. (where=(d_stfid ne o_stfid));
        title "imputation rate of tran for &met.";
        tables predict_tran;
    run;
    proc means data=OUTDATA.proxmeas_&met. 
        (where=((predict_tran=0) and (d_stfid ne o_stfid)));
        title "summary of reported tran travel time for &met.";
        var dist tran speed_tran
        %if ("&use_prox_time_mpo_auto."="1") %then %do;
        auto speed_auto 
        better_tran best speed_best;
        %end;
        ;
    run;
    %end;          
    
    * end evaluation of mpo proximity;
    %end;

    * alternate long term fix for missing travel time data (not done);
    * identify each missing time, make a dataset;
    * make a dataset sorted by work geo and distance, retain shortest distance obs;
    * merge shortest distance obs with travel time data by home work;
    * add distance to donor tract to travel time;
    * merge back on to primary dataset;

    * merge CTPP travel time to MPO time;
    %if ("&use_ctpp_time."="1") %then %do;
    * merge on CTPP variables;
    libname DATACTPP "&dirdata./input_data/ctpp_pu/&st." access=readonly;
    
    * first, read in CTPP for state;
    * TABLE306 - Means of Transportation (18);
    * tab306x02 Drove alone;
    * tab306x08 Bus or trolley bus;
    * tab306x09 Streetcar or trolley car;
    * tab306x10 Subway or elevated;
    * tab306x11 Railroad;
    * tab306x12 Ferryboat;
    * TABLE309 - Median travel time by Means of Transportation (8) and Time Leaving Home (4);
    * tab309x06: Drove alone, Time leaving home 5:00 a.m. to 8:59 a.m.;
    * tab309x18: Bus or trolley bus, Time leaving home 5:00 a.m. to 8:59 a.m.;
    * tab309x22: Streetcar, subway, rail or ferry, Time leaving home 5:00 a.m. to 8:59 a.m.;
    data ctpp_time;
        length o_stfid d_stfid $11;
        set DATACTPP.&st._sumlv140 (keep=state3 county tract qpowst qpowco qpowtract
                                    tab306x02 tab306x08 tab306x09 tab306x10 tab306x11
                                    tab309x06 tab309x18 tab309x22);
        * geography;
        o_stfid=substr(state3,2,2)||county||tract;
        d_stfid=substr(qpowst,2,2)||qpowco||qpowtract;
        drop state3 county tract qpowst qpowco qpowtract;
        * means of transport;
        ctpp_auto_n=tab306x02;
        ctpp_tran_n=tab306x08+tab306x09+tab306x10+tab306x11;
        ctpp_tran_n_bus=tab306x08;
        ctpp_tran_n_rail=tab306x09+tab306x10+tab306x11;
        drop tab306x02 tab306x08 tab306x09 tab306x10 tab306x11;
        * calculate travel time;
        ctpp_auto=.;
        if (ctpp_auto_n>0 and tab309x06>0) then ctpp_auto=tab309x06;
        ctpp_tran=.;
        if ((tab309x18>0 and ctpp_tran_n_bus>0) and (tab309x22>0 and ctpp_tran_n_rail>0)) 
        then ctpp_tran=(tab309x18*ctpp_tran_n_bus+tab309x22*ctpp_tran_n_rail)/ctpp_tran_n;
        if ((tab309x18>0 and ctpp_tran_n_bus>0) and (tab309x22=0 or ctpp_tran_n_rail=0)) 
        then ctpp_tran=tab309x18;
        if ((tab309x18=0 or ctpp_tran_n_bus=0) and (tab309x22>0 and ctpp_tran_n_rail>0)) 
        then ctpp_tran=tab309x22;
        drop tab309x06 tab309x18 tab309x22;
        * only retain records with at least some flows;
        if (ctpp_auto_n>0 or ctpp_tran_n>0);
    run;
    
    * merge to mpo data;
    proc sort data=ctpp_time;
        by d_stfid o_stfid;
    run;
    data OUTDATA.proxmeas_&met.;
        merge OUTDATA.proxmeas_&met (in=in_prox)
              ctpp_time (in=in_ctpp);
        by d_stfid o_stfid;
        * keep origins and destinations in proximity file, in analysis area;
        if (in_prox);
        * indicator for match to ctpp;
        if (in_ctpp) then do;
            flow_ctpp=1;
            * indicator of sufficient auto data;
            if ((ctpp_auto_n ge 10) and (ctpp_auto gt 0)) then flow_ctpp_auto=1;
            else flow_ctpp_auto=0;
            * indicator of sufficient transit data;
            if ((ctpp_tran_n ge 10) and (ctpp_tran gt 0)) then flow_ctpp_tran=1;
            else flow_ctpp_tran=0;
        end;
        else flow_ctpp=0;
    run;
                                    
    * analysis of ctpp match;
    proc freq data=OUTDATA.proxmeas_&met.;
        title "match rate to ctpp for &met.";
        table flow_ctpp;
    run;

    * analysis of ctpp match - auto;
    proc freq data=OUTDATA.proxmeas_&met.;
        title "auto match rate to ctpp and prediction rate for &met.";
        table flow_ctpp_auto*predict_auto;
    run;
    * analysis of ctpp match - transit;
    proc freq data=OUTDATA.proxmeas_&met.;
        title "tran match rate to ctpp and prediction rate for &met.";
        table flow_ctpp_tran*predict_tran;
    run;
    
    * means comparison - auto;
    proc means data=OUTDATA.proxmeas_&met. (where=(flow_ctpp_auto=1 and predict_auto=0));
        title "summary of ctpp v mpo auto time for &met.";
        var ctpp_auto auto;
    run;
    * means comparison - transit;
    proc means data=OUTDATA.proxmeas_&met. (where=(flow_ctpp_tran=1 and predict_tran=0));
        title "summary of ctpp v mpo tran time for &met.";
        var ctpp_tran tran;
    run;

    /* figure - auto;
    filename outgraph "&dirprogvint./loglst/mpo_ctpp_auto_&met..gif";
    goptions gsfname=outgraph dev=gif373;
    title "plot of auto times, ctpp v mpo for &met.";
    proc gplot data=OUTDATA.proxmeas_&met. (where=(flow_ctpp_auto=1));
        plot auto*ctpp_auto / haxis=axismpo vaxis=axisctpp;
        * define axis characteristics;
        axismpo  label=(angle=90 "MPO time");
        axisctpp label=("CTPP time");
    run;
    quit;
    * figure - transit;
    filename outgraph "&dirprogvint./loglst/mpo_ctpp_tran_&met..gif";
    goptions gsfname=outgraph dev=gif373;
    title "plot of transit times, ctpp v mpo for &met.";
    proc gplot data=OUTDATA.proxmeas_&met. (where=(flow_ctpp_tran=1));
        plot tran*ctpp_tran  / haxis=axismpo vaxis=axisctpp;
        * define axis characteristics;
        axismpo  label=(angle=90 "MPO time");
        axisctpp label=("CTPP time");
    run;
    quit;*/
        
    * estimation of ctpp travel time based on mpo times - auto;
    proc reg data=OUTDATA.proxmeas_&met. (where=(flow_ctpp_auto=1 and predict_auto=0));
        title "regression of auto, ctpp travel time for &met.";
        model ctpp_auto=auto;
    run;
    * estimation of ctpp travel time based on mpo times - transit;
    proc reg data=OUTDATA.proxmeas_&met. (where=(flow_ctpp_tran=1 and predict_tran=0));
        title "regression of tran, ctpp travel time for &met.";
        model ctpp_tran=tran;
    run;
            
    * end merge CTPP travel time to MPO time;
    %end;

    
%end; * end proximity;

/* create frame of attributes to use */
/* run OTM attributes --------------------------------------------------------------------- */
%if (&run_otm_attrib.=1) %then %do;
    %macro otmdatapull(otmarea,otmpref,otmyear,otmvars,outname);
        * identify tracts used for city wac_what_b_dom_2005.sas7bdat;
        data trt_attrib (drop=&otmpref._stfid);
            length stfid $11;
            * read in multiple earnings segments here and combine;
            set OTM.&otmarea._what_b_dom_&otmyear. 
                (keep=&otmpref._stfid &otmvars. 
                   where=(owner='B3' 
                   and (input(substr(&otmpref._stfid,1,2),2.) = &stcode.)
                   and (substr(&otmpref._stfid,3,3) in (&cty.))
                   ));
            stfid=substr(&otmpref._stfid,1,11);
            * year;
            year=&otmyear.;
            * for most cases, limit to low to moderate earning jobs;
            if (earn_c='1' or earn_c='2') then do;
                  * rename or reformat otm variables;
                  all=jobs;
                  * Goods 11, 21, 22, 23, 31-33 Distribution 42, 48-49;
                  if (naics in ('01','02','03','04','05','06','08')) then sc1=jobs;
                  else sc1=0;
                  * Retail 44-45 Service to customers 71, 72 Other services 56, 81;
                  if (naics in ('07','17','18','14','19')) then sc2=jobs;
                  else sc2=0;
                  * Professional 51, 54, 55 FIRE 52, 53;
                  if (naics in ('09','12','13','10','11')) then sc3=jobs;
                  else sc3=0;
                  * Public 61, 92;
                  if (naics in ('15','20')) then sc4=jobs;
                  else sc4=0;
                  * Health 62;
                  if (naics in ('16')) then sc5=jobs;
                  else sc5=0;
                  * White race/ethnicity;
                  if (race='1' and ethnicity='N') then gp1=jobs;
                  else gp1=0;
                  * Black race/ethnicity;
                  if (race='2' and ethnicity='N') then gp2=jobs;
                  else gp2=0;
                  * Hispanic race/ethnicity;
                  if (ethnicity='H') then gp3=jobs;
                  else gp3=0;
                  * Other race/ethnicity;
                  if (ethnicity='N' and race>'2') then gp4=jobs;
                  else gp4=0;
            end;
            * all earnings levels;
            eall=jobs;
        run;
        proc summary data=trt_attrib nway;
            class stfid year;
            var &use_otmvars2.;
            output out=OUTDATA.attrib_otm_&outname._&met._&otmyear. (drop=_TYPE_ _FREQ_) sum=;
            * use / autoname if anything other than sum;
        run;
        proc contents data=OUTDATA.attrib_otm_&outname._&met._&otmyear.;
            title "check variables from otm &outname._&met._&otmyear.";
        run;
    %mend otmdatapull;
    * load in OTM rac and wac and output by tract for cities we care about;
    * Note that these are summaries of the WHATB, not directly from LODES/OTM;
    * work;
    %if (&use_otmjob.=1) %then %do;
        %otmdatapull(wac,w,&otmyear.,&use_otmvars.,job)
    %end;
    * home;
    %if (&use_otmlab.=1) %then %do;
        %otmdatapull(rac,h,&otmyear.,&use_otmvars.,lab)
    %end;
%end; * end attributes;


/* run accessibility ---------------------------------------------------------------------- */
%if (&run_accessibility.=1) %then %do;
/* merge frame of attributes to tracts and calculate measures */
    %macro merge_attrib_otm(geomerge,att);
        %let otmvarsrename=;
        %let otmvar=1;
        %do %until("%scan(&use_otmvars2.,&otmvar.)"="");
            %let otmvarname=%scan(&use_otmvars2.,&otmvar.);
            %let otmvarsrename=&otmvarsrename. &otmvarname.=&geomerge._&att._&otmvarname.;
            %let otmvar=%eval(&otmvar+1);
        %end;
        OUTDATA.attrib_otm_&att._&met._&otmyear. (in=in_&geomerge. rename=(stfid=&geomerge._stfid &otmvarsrename.))
    %mend merge_attrib_otm;
    * merge attributes by workplace or destination;
    data trt_merge_attrib;
        merge OUTDATA.proxmeas_&met. (in=in_prox)
              %if (&use_otmjob.=1) %then %do;
                  %merge_attrib_otm(d,job)
              %end;
              %if (&use_otmlab.=1) %then %do;
                  %merge_attrib_otm(d,lab)
              %end;
              ;
        by d_stfid;
        if in_prox;       
    run;
    * sort by home then work before sum by home (could also sort by theta here);
    proc sort data=trt_merge_attrib;
        by o_stfid d_stfid;
    run;
    
    * export for further analysis (d_job_c000);
    %if (&export_access.=1) %then %do;
    data OUTDATA.trt_merge_attrib_&met._&otmyear. 
        (keep=o_stfid d_stfid dist auto /*tran d_job_ce01*/);
        set trt_merge_attrib;
    run;
    proc export data=OUTDATA.trt_merge_attrib_&met._&otmyear.
        outfile="&dirdatavint./trt_merge_attrib_&met._&otmyear..csv"
        dbms=csv
        replace;
    run;
    %end;

    * merge attributes by residence or origin;
    data trt_merge_attrib;
        merge trt_merge_attrib (in=in_prox)
              /*%if (&use_otmjob.=1) %then %do;
                  %merge_attrib_otm(o,job)
              %end;*/
              %if (&use_otmlab.=1) %then %do;
                  %merge_attrib_otm(o,lab)
              %end;
              ;
        by o_stfid;
        if in_prox;       
    run;
    
    * macros specify a calculation, proximity measure, attribute, and modifier;
    * macro to use in creating new variable names to be printed, for length statement;
    %macro print_varlist(ren_calc,ren_prox,ren_att,ren_mod);
        %let otmvars_newname=;
        %let otmvar=1;
        %do %until("%scan(&use_otmvars2.,&otmvar.)"="");
            %let otmvarname=%scan(&use_otmvars2.,&otmvar.);
            %let otmvars_newname=&otmvars_newname. &ren_calc._&ren_prox._&ren_att._&ren_mod._&otmvarname.;
            %let otmvar=%eval(&otmvar+1);
        %end;
        &otmvars_newname.
    %mend print_varlist;
    * macro to use in creating new variable names use in another macro;
    %macro print_varlist_attrib(ren_att);
        %let otmvars_newname=;
        %let otmvar=1;
        %do %until("%scan(&use_otmvars2.,&otmvar.)"="");
            %let otmvarname=%scan(&use_otmvars2.,&otmvar.);
            %let otmvars_newname=&otmvars_newname. &ren_att._&otmvarname.;
            %let otmvar=%eval(&otmvar+1);
        %end;
        &otmvars_newname.
    %mend print_varlist_attrib;
    * array statement for new variables;
    %macro array_setup(calc,prox,att,mod);
        array array_&calc._&prox._&att._&mod. (&use_otmvars2_count.) %print_varlist(&calc.,&prox.,&att.,&mod.);
    %mend array_setup;
    * index calculation statement;
    %macro spec_index(prox,att,mod);
        if (&prox.=.) then &prox.=0;
        do i=1 to &use_otmvars2_count.;
            if (&prox.=. or array_&att.(i)=.) then array_index_&prox._&att._&mod.(i)=0;
            * modify so that variable does not include decimal point;
            else array_index_&prox._&att._&mod.(i)=array_&att.(i)/((1+&prox.)**(&mod.*0.01));
        end;
    %mend spec_index;
    * index modified calculation statement;
    %macro spec_indmd(prox,att,mod);
        if (&prox.=.) then &prox.=0;
        do i=1 to &use_otmvars2_count.;
            if (&prox.=. or array_&att.(i)=.) then array_indmd_&prox._&att._&mod.(i)=0;
            * modify so that variable does not include decimal point;
            else if (&prox.<&mod.) then array_indmd_&prox._&att._&mod.(i)=array_&att.(i);
            * after a certain range, begin discounting by surplus range;
            else if (&prox.<60) then
            array_indmd_&prox._&att._&mod.(i)=array_&att.(i)/((&prox.-(&mod.-1))**(1));
            else array_indmd_&prox._&att._&mod.(i)=0;
        end;
    %mend spec_indmd;
    * index modified calculation statement, exponential;
    %macro spec_indxx(prox,att,mod,exp);
        if (&prox.=.) then &prox.=0;
        do i=1 to &use_otmvars2_count.;
            if (&prox.=. or array_&att.(i)=.) then array_indxx_&prox._&att._&mod.(i)=0;
            * modify so that variable does not include decimal point;
            else if (&prox.<&mod.) then array_indxx_&prox._&att._&mod.(i)=array_&att.(i);
            * after a certain range, begin discounting by surplus range;
            else if (&prox.<60) then
            array_indxx_&prox._&att._&mod.(i)=array_&att.(i)/(exp(&exp.*(max(1,&prox.-&mod.))));
            else array_indxx_&prox._&att._&mod.(i)=0;
        end;
    %mend spec_indxx;
    * index calculation statement;
    %macro spec_index_poly(prox,att,mod);
        if (&prox.=.) then &prox.=0;
        do i=1 to &use_otmvars2_count.;
            if (&prox.=. or array_&att.(i)=.) then array_index_&prox._&att._&mod.(i)=0;
            * modifier is just a stand in here, divide by exp of constant to normalize;
            else array_index_&prox._&att._&mod.(i)=array_&att.(i)*exp(3-0.125*&prox.+0.001*(&prox.**2))*(1/exp(3));
        end;
    %mend spec_index_poly;
    * rings calculation statement;
    %macro spec_rings(prox,att,mod);
        if (&prox. ge &mod.) then do;
            do i=1 to &use_otmvars2_count.;
                array_rings_&prox._&att._&mod.(i)=0;
            end;
        end;
        if (&prox. lt &mod.) then do;
            do i=1 to &use_otmvars2_count.;
                if (&prox.=. or array_&att.(i)=.) then array_rings_&prox._&att._&mod.(i)=0;
                else array_rings_&prox._&att._&mod.(i)=array_&att.(i);
            end;
        end;
    %mend spec_rings;
    * donut calculation statement;
    %macro spec_donut(prox,att,mod,low);
        if (&prox. ge &mod.) then do;
            do i=1 to &use_otmvars2_count.;
                array_donut_&prox._&att._&mod.(i)=0;
            end;
        end;
        if (&prox. ge &low. and &prox. lt &mod.) then do;
            do i=1 to &use_otmvars2_count.;
                if (&prox.=. or array_&att.(i)=.) then array_donut_&prox._&att._&mod.(i)=0;
                else array_donut_&prox._&att._&mod.(i)=array_&att.(i);
            end;
        end;
    %mend spec_donut;
    * create measures - index and rings - in one data step;
    data access_calc; * (drop=job_: lab_:);
        * define new variables, list combinations most interested in,;
        * must be consistent in arrays and body below;
        length
            %print_varlist(indxx,auto,d_job,15)
            %print_varlist(indx,auto,o_lab,15)
            %print_varlist(indxx,auto,d_job,10)
            %print_varlist(indxx,auto,o_lab,10)
            %print_varlist(indxx,auto,d_job,5)
            %print_varlist(indxx,auto,o_lab,5)
            %print_varlist(indxx,tran,d_job,10)
            %print_varlist(indxx,tran,o_lab,10)
            %print_varlist(indxx,dist,d_job,1)
            %print_varlist(indxx,dist,o_lab,1)
            %print_varlist(rings,dist,d_job,3)
            %print_varlist(rings,dist,o_lab,3)
            5.;
        * dataset of proximity measures and attributes;
        set trt_merge_attrib;
        * define year of otm data to be integrated (even though year by year processing);
        year=&otmyear.;
        * define arrays of input attributes;
        %if (&use_otmjob.=1) %then %do;
        array array_d_job(&use_otmvars2_count.) %print_varlist_attrib(d_job);
        %end;
        %if (&use_otmlab.=1) %then %do;
        array array_d_lab(&use_otmvars2_count.) %print_varlist_attrib(d_lab);
        array array_o_lab(&use_otmvars2_count.) %print_varlist_attrib(o_lab);
        %end;
        * define arrays of new variables;
        %array_setup(indxx,auto,d_job,15);
        %array_setup(indxx,auto,o_lab,15);
        %array_setup(indxx,auto,d_job,10);
        %array_setup(indxx,auto,o_lab,10);
        %array_setup(indxx,auto,d_job,5);
        %array_setup(indxx,auto,o_lab,5);
        %array_setup(indxx,tran,d_job,10);
        %array_setup(indxx,tran,o_lab,10);
        %array_setup(indxx,dist,d_job,1);
        %array_setup(indxx,dist,o_lab,1);
        %array_setup(rings,dist,d_job,3);
        %array_setup(rings,dist,o_lab,3);
        * specify what to do if a missing value for a time;
        * make far away, could be on input side;
        * divide attribute by proximity normalizer;
        %spec_indxx(auto,d_job,15,0.1);
        %spec_indxx(auto,o_lab,15,0.1);
        %spec_indxx(auto,d_job,10,0.1);
        %spec_indxx(auto,o_lab,10,0.1);
        %spec_indxx(auto,d_job,5,0.1);
        %spec_indxx(auto,o_lab,5,0.1);
        %spec_indxx(tran,d_job,10,0.1);
        %spec_indxx(tran,o_lab,10,0.1);
        %spec_indxx(dist,d_job,1,1);
        %spec_indxx(dist,o_lab,1,1);
        %spec_rings(dist,d_job,3);
        %spec_rings(dist,o_lab,3);
        drop i;
    run;
    * output small file for analysis;
    data OUTDATA.access_calc;
        set access_calc (obs=1000);
    run;

    * these variable are now called rings or index, but have not yet been summed up;
    * do initial sum of competing searchers able to reach a destination, used below;
    proc sort data=access_calc;
        by d_stfid;
    run;
    proc summary data=access_calc nway;
        by d_stfid year;
        var indxx_auto_o_lab: indxx_dist_o_lab: rings_dist_o_lab:
            ;
        output out=OUTDATA.compete_&met._&otmyear. (drop=_TYPE_ _FREQ_) sum=;
    run;
    * multiply jobs at workplace by labor force that can commute there; 
    * after the next proc sum, divide these products by total jobs accessible from home;
    * the result is a job weighted average of labor force competition for those jobs;
    %macro ringscalc(prox,ringval,ringvar);
        rings_&prox._o_labw_&ringval._&ringvar.=
        max(0,rings_&prox._o_lab_&ringval._&ringvar.)*max(0,rings_&prox._d_job_&ringval._&ringvar.);
    %mend ringscalc;
    %macro donutcalc(prox,ringval,ringvar);
        donut_&prox._o_labw_&ringval._&ringvar.=
        max(0,rings_&prox._o_lab_&ringval._&ringvar.)*max(0,donut_&prox._d_job_&ringval._&ringvar.);
    %mend donutcalc;
    %macro indmdcalc(prox,indmdval,indmdvar);
        indmd_&prox._o_labw_&indmdval._&indmdvar.=
        max(0,indmd_&prox._o_lab_&indmdval._&indmdvar.)*max(0,indmd_&prox._d_job_&indmdval._&indmdvar.);
    %mend indmdcalc;
    %macro indxxcalc(prox,indxxval,indxxvar);
        indxx_&prox._o_labw_&indxxval._&indxxvar.=
        max(0,indxx_&prox._o_lab_&indxxval._&indxxvar.)*max(0,indxx_&prox._d_job_&indxxval._&indxxvar.);
    %mend indxxcalc;
    %macro indxxcalc2(prox,indxxval,indxxvar,prox2);
        indxx_&prox._o_labw_&indxxval._&indxxvar.=
        max(0,indxx_&prox2._o_lab_&indxxval._&indxxvar.)*max(0,indxx_&prox._d_job_&indxxval._&indxxvar.);
    %mend indxxcalc2;
    * merge back with other data (drop original, non-summed searchers;
    data access_calc_&met._pre;
        merge access_calc (in=workplaces 
                           drop=indxx_auto_o_lab: indxx_dist_o_lab: rings_dist_o_lab:)
              OUTDATA.compete_&met._&otmyear.;
        by d_stfid;
        * only do calculations within each workplace;
        if (workplaces);
        * first ring is just a circle;
        %let otmvar=1;
        %do %until("%scan(&use_otmvars2.,&otmvar.)"="");
            %let otmvarname=%scan(&use_otmvars2.,&otmvar.);
            %indxxcalc(auto,15,&otmvarname.)
            %indxxcalc(auto,10,&otmvarname.)
            %indxxcalc(auto,5,&otmvarname.)
            %indxxcalc2(tran,10,&otmvarname.,auto)
            %indxxcalc(dist,1,&otmvarname.)
            %ringscalc(dist,3,&otmvarname.)
            %let otmvar=%eval(&otmvar+1);
        %end;
    run;
        
    * sum up by home tract;
    proc sort data=access_calc_&met._pre;
        by o_stfid d_stfid;
    run;
    proc summary data=access_calc_&met._pre nway;
        by o_stfid year;
        var indxx_: rings_:;
        output out=access_&met._intwrk (drop=_TYPE_ _FREQ_) sum=;
    run;

    * macro for calculating weighted average;
    * (in case of no jobs, use 1 here to prevent missing values, will be 0/1=0);
    %macro weightavg(calc,prox,att,mod,otmvar);
        &calc._&prox._&att._&mod._&otmvar.=
        &calc._&prox._&att.w_&mod._&otmvar./max(1,&calc._&prox._d_job_&mod._&otmvar.);
        drop &calc._&prox._&att.w_&mod._&otmvar.;
    %mend weightavg;
    * macro for normalizing job labor ratio (redone in Stata code);
    %macro rationorm(calc,prox,att,mod,otmvar);
        &calc._&prox._ratio_&mod._&otmvar.=
        (max(1,&calc._&prox._d_job_&mod._&otmvar.)-max(1,&calc._&prox._&att._&mod._&otmvar.))/
        (0.5*(max(1,&calc._&prox._d_job_&mod._&otmvar.)+max(1,&calc._&prox._&att._&mod._&otmvar.)));
        if (&calc._&prox._ratio_&mod._&otmvar.=.) then &calc._&prox._ratio_&mod._&otmvar.=0;
    %mend rationorm;
    * final computations;
    data OUTDATA.access_&met._&otmyear.;
        retain o_stfid year;
        set access_&met._intwrk;
        * job weighted average of competing searchers or labor force;
        * for use as a separate variable of competing searchers, or input to ratios;
        %macro weightavgvars(avgvar);
            %weightavg(indxx,auto,o_lab,15,&avgvar.);
            %weightavg(indxx,auto,o_lab,10,&avgvar.);
            %weightavg(indxx,auto,o_lab,5,&avgvar.);
            %weightavg(indxx,tran,o_lab,10,&avgvar.);
            %weightavg(indxx,dist,o_lab,1,&avgvar.);
            %weightavg(rings,dist,o_lab,3,&avgvar.);
        %mend weightavgvars;
        %weightavgvars(all)
        %weightavgvars(sc1)
        %weightavgvars(sc2)
        %weightavgvars(sc3)
        %weightavgvars(sc4)
        %weightavgvars(sc5)
        %weightavgvars(gp1)
        %weightavgvars(gp2)
        %weightavgvars(gp3)
        %weightavgvars(gp4)
        %weightavgvars(eall)
        * job to labor force ratio, normalized ratio, 0 to 2;
        %macro rationormvars(avgvar);
            %rationorm(indxx,auto,o_lab,15,&avgvar.);
            %rationorm(indxx,auto,o_lab,10,&avgvar.);
            %rationorm(indxx,auto,o_lab,5,&avgvar.);
            %rationorm(indxx,tran,o_lab,10,&avgvar.);
            %rationorm(indxx,dist,o_lab,1,&avgvar.);
            %rationorm(rings,dist,o_lab,3,&avgvar.);
        %mend rationormvars;
        %rationormvars(all)
        %rationormvars(sc1)
        %rationormvars(sc2)
        %rationormvars(sc3)
        %rationormvars(sc4)
        %rationormvars(sc5)
        %rationormvars(gp1)
        %rationormvars(gp2)
        %rationormvars(gp3)
        %rationormvars(gp4)
        %rationormvars(eall)
    run;
                                        
%end; * end accessibility;


/* run summary ---------------------------------------------------------------------- */
%if (&run_summary.=1) %then %do;

    proc means data=OUTDATA.attrib_otm_job_&met._&otmyear. sum mean n nway;
        title "tract level job summary, for &met.";    
        var all
            ;
        output out=job_means sum=sum_job_all mean=mean_job_all n=n_job_all;
    run;
    proc means data=OUTDATA.attrib_otm_lab_&met._&otmyear. sum mean n nway;
        title "tract level labor summary, for &met.";    
        var all
            ;
        output out=job_means sum=sum_lab_all mean=mean_lab_all n=n_lab_all;        
    run;

    proc means data=OUTDATA.access_&met._&otmyear. n mean min max nway;
        class year;
        title "tract and metro level spatial mismatch summary, for &met., by year";
        var indxx_auto_d_job_10_all
            indxx_auto_o_lab_10_all
            indxx_auto_ratio_10_all
            ;
        output out=access_means n= mean= min= max= /autoname;        
    run;
            
    proc means data=OUTDATA.access_&met._&otmyear. n mean min max nway;
        title "tract and metro level spatial mismatch summary, for &met.";
        var indxx_auto_d_job_15_all
            indxx_auto_o_lab_15_all
            indxx_auto_d_job_10_all
            indxx_auto_o_lab_10_all
            indxx_tran_d_job_10_all
            indxx_tran_o_lab_10_all
            indxx_dist_d_job_1_all
            indxx_dist_o_lab_1_all
            rings_dist_d_job_3_all
            rings_dist_o_lab_3_all
            indxx_auto_ratio_15_all
            indxx_auto_ratio_10_all
            indxx_tran_ratio_10_all
            indxx_dist_ratio_1_all
            ;
        output out=access_means n= mean= min= max= /autoname;        
    run;

%end;


%mend tract_access; *end of macro run in primary code;

                
/* run census 2000 SF3 ---------------------------------------------------------------------- */
%macro cen2000_sf3;

    * home census 2000;
    data OUTDATA.census2000_sf3 (keep=stfid pop pop25 pop25nco pop16 pop16lab);
        length stfid $11 ss $2 ccc $3 tttttt $6 pop pop25 pop25nco pop16 pop16lab 5.;
        set CEN2000.census2000_sf3;
        ss=put(state,z2.);
        ccc=put(county,z3.);
        tttttt=put(tract,z6.);
        stfid=ss||ccc||tttttt;
        pop=max(0,p001001);
        pop25=max(0,p037001);
        pop25nco=pop25-sum(P037015+P037016+P037017+P037018)-sum(P037032+P037033+P037034+P037035);
        pop16=max(0,P043001);
        pop16lab=sum(P043003+P043010);
        label
            pop="Census 2000 tract population"
            pop25="Census 2000 tract population age 25 and up"
            pop25nco="Census 2000 tract population age 25 and up, no college"
            pop16="Census 2000 tract population age 16 and up"
            pop16lab="Census 2000 tract population age 16 and up, in labor force"
        ;
    run;
    proc means data=OUTDATA.census2000_sf3;
        title "SF3 Census 2000";    
        var pop pop25 pop25nco pop16 pop16lab;
    run;

    * merge attributes by residence or origin;
    * index modified calculation statement, exponential;
    %macro spec_indxx_sf3(prox,att,mod,exp,popvar);
        if (&prox.=.) then &prox.=0;
        if (&prox.=. or &popvar.=.) then indxx_&prox._&att._&mod._&popvar.=0;
        * modify so that variable does not include decimal point;
        else if (&prox.<&mod.) then indxx_&prox._&att._&mod._&popvar.=&popvar.;
        * after a certain range, begin discounting by surplus range;
        else if (&prox.<60) then
        indxx_&prox._&att._&mod._&popvar.=&popvar./(exp(&exp.*(max(1,&prox.-&mod.))));
        else indxx_&prox._&att._&mod._&popvar.=0;
    %mend spec_indxx_sf3;
    * assemble proximity data (limit time to 40 min or something);
    data proxmeas;
        set
        %let metord=1;
        %do %until("%scan(&metrolist.,&metord.)"="");
            %let met=%scan(&metrolist.,&metord.);
            %let stmet=%scan(&stmetrolist.,&metord.);
            OUTDATA.proxmeas_&met.
            (keep=d_stfid o_stfid auto tran dist
            where=((auto lt 60) or (tran lt 60) or (dist lt 60)))
            %let metord=%eval(&metord+1);
        %end;
        ;
    run;
    proc sort data=proxmeas;
        by o_stfid d_stfid;
    run;
    data proxmeas;
        merge proxmeas (in=in_prox)
              OUTDATA.census2000_sf3 (rename=(stfid=o_stfid))
              ;
        by o_stfid;
        if in_prox;
        %spec_indxx_sf3(auto,o_lab,10,0.1,pop16lab);
        %spec_indxx_sf3(auto,o_lab,10,0.1,pop25nco);
    run;
    proc sort data=proxmeas;
        by d_stfid o_stfid;
    run;

    * output small file for analysis;
    *data OUTDATA.sf3_calc_check;
    *    set proxmeas (obs=1000);
    *run;

    * do initial sum of competing searchers able to reach a destination, used below;
    proc summary data=proxmeas nway;
        by d_stfid;
        var indxx_auto_o_lab:;
        output out=OUTDATA.compete_cen2000 (drop=_TYPE_ _FREQ_) sum=;
    run;
    proc means data=OUTDATA.compete_cen2000;
        title "SF3 Census 2000, access";    
        var indxx_auto_o_lab:;
    run;
    
%mend cen2000_sf3;     
                
                
/* run census 2000 SF3 summary ------------------------------------------------------------ */
%macro cen2000_sf3_sum;
    
    * aggregate job access;
    * initialize code macros that will be used below;
    %let metord = 1;
    * cycle through each city, choosing the appropriate state, state codes, and counties;
    %do %until("%scan(&metrolist.,&metord.)"="");
        %let met=%scan(&metrolist.,&metord.);
        %let st=%scan(&stlist.,&metord.);
        %let stmet=%scan(&stmetrolist.,&metord.);
        %let stcode=%scan(&stcodelist.,&metord.);
        %let cty=%scan(%quote(&countylist.),&metord.,V);
        %put  &met. &st. &cty.;

        proc append base=access_met
                    data=OUTDATA.access_&met. 
                        (keep=o_stfid year indxx_auto_d_job_10_all indxx_auto_o_lab_10_all);
        run;

        %let metord=%eval(&metord+1);
    %end;            
                
    * summarize;
    proc sort data=OUTDATA.census2000_sf3 (keep=stfid pop pop25 pop25nco pop16 pop16lab 
              rename=(stfid=h_stfid)) out=cen;
        by h_stfid;
    run;
    proc sort data=access_met (keep=o_stfid year indxx_:
              rename=(o_stfid=h_stfid)) out=acc;
        by h_stfid;
    run;
    data acc;
        merge acc (in=in_acc)
              cen (in=in_cen)
              ;
        by h_stfid;
        if (in_acc);
        ratio_auto_10_all=(max(1,indxx_auto_d_job_10_all)-max(1,indxx_auto_o_lab_10_all))/
                (0.5*(max(1,indxx_auto_d_job_10_all)+max(1,indxx_auto_o_lab_10_all)));
    run;
    proc means data=acc n p10 median p90;
        title "JA population weighted";
        var ratio_auto_10_all;
        weight pop;
    run;
    proc means data=acc n p10 median p90;
        title "JA labor market weighted";
        var ratio_auto_10_all;
        weight pop16lab;
    run;
    proc means data=acc n p10 median p90;
        title "JA no college weighted";
        var ratio_auto_10_all;
        weight pop25nco;
    run;
    proc means data=acc n p10 median p90;
        class year;
        title "JA population weighted";
        var ratio_auto_10_all;
        weight pop;
    run;
    proc means data=acc n p10 median p90;
        class year;
        title "JA labor market weighted";
        var ratio_auto_10_all;
        weight pop16lab;
    run;
    proc means data=acc n p10 median p90;
        class year;
        title "JA no college weighted";
        var ratio_auto_10_all;
        weight pop25nco;
    run;
    
%mend cen2000_sf3_sum;


* for all states together;
%macro us_access;

/* run super-summary ---------------------------------------------------------------------- */
%if (&run_supersummary.=1) %then %do;

    * initialize code macros that will be used below;
    %let metord = 1;
    * cycle through each city, choosing the appropriate state, state codes, and counties;
    %do %until("%scan(&metrolist.,&metord.)"="");
        %let met=%scan(&metrolist.,&metord.);
        %let st=%scan(&stlist.,&metord.);
        %let stmet=%scan(&stmetrolist.,&metord.);
        %let stcode=%scan(&stcodelist.,&metord.);
        %let cty=%scan(%quote(&countylist.),&metord.,V);
        %put  &met. &st. &cty.;

        * grab next metro file;
        data proxmeas;
            set OUTDATA.proxmeas_&met. 
            (keep=o_stfid d_stfid dist auto tran 
            predict_auto speed_auto predict_tran speed_tran
            flow_ctpp_auto flow_ctpp_tran ctpp_auto ctpp_tran
            );
        run;

        * append to national file;
        proc append base=proxmeas_national data=proxmeas;
        run;

        %let metord=%eval(&metord+1);
    %end;

    * resort by stfid order of home;        
    proc sort data=proxmeas_national out=OUTDATA.proxmeas_national;
        by o_stfid d_stfid;
    run;
            
    * compress LODES od file to tract level;
    data flows_od (keep=o_stfid d_stfid jobs);
        set OTMFLOW.od_us_2005_jt03_s000_1 (where=(substr(h_geocode,1,2) in (&stcodelistq.)));
        o_stfid=substr(h_geocode,1,11);
        d_stfid=substr(w_geocode,1,11);
    run; 
    proc summary data=flows_od nway;
        class o_stfid d_stfid;
        var jobs;
        output out=flows_od (drop=_TYPE_ _FREQ_) sum=;
    run;
        
    * merge to attribute data;
    data OUTDATA.proxmeas_national;
        length o_stfid d_stfid $11;
        merge OUTDATA.proxmeas_national (in=in_prox)
              flows_od (in=in_flow)
              ;
        by o_stfid d_stfid;
        * keep only if in prox file;
        if (in_prox);
        * set jobs to zero if no jobs;
        if (in_flow ne 1) then jobs=0;
        * calculate ctpp speeds;
        if (flow_ctpp_auto=1) then do;
            if (o_stfid=d_stfid) then speed_ctpp_auto=.;
            else speed_ctpp_auto=(dist/ctpp_auto)*60;
        end;
        if (flow_ctpp_tran=1) then do;
            if (o_stfid=d_stfid) then speed_ctpp_tran=.;
            else speed_ctpp_tran=(dist/ctpp_tran)*60;
        end;
    run;    
            
    * summarize auto non-impute; 
    proc means data=OUTDATA.proxmeas_national 
        (where=(o_stfid ne d_stfid and predict_auto=0)) n median mean std min max;
        title "summarize auto non-impute";
        weight jobs;
        var dist auto speed_auto
            ;
        output out=access_means_autompo n= median= mean= std= min= max= /autoname;        
    run;
        
    * summarize transit non-impute; 
    proc means data=OUTDATA.proxmeas_national 
        (where=(o_stfid ne d_stfid and predict_tran=0 and speed_tran gt 3)) n median mean std min max;
        title "summarize transit has transit";
        weight jobs;
        var dist tran speed_tran
            ;
        output out=access_means_tranmpo n= median= mean= std= min= max= /autoname;        
    run;
        
    * summarize transit non-impute; 
    proc means data=OUTDATA.proxmeas_national 
        (where=(o_stfid ne d_stfid and speed_tran eq 3)) n median mean std min max;
        title "summarize transit no transit";
        weight jobs;
        var dist tran speed_tran
            ;
        output out=access_means_tranimp n= median= mean= std= min= max= /autoname;        
    run;
        
    * summarize all tracts with flows; 
    proc means data=OUTDATA.proxmeas_national 
        (where=(o_stfid ne d_stfid)) n median mean std min max;
        title "summarize all tracts with flows";
        weight jobs;
        var dist auto tran speed_auto speed_tran
            ;
        output out=access_means_proxall n= median= mean= std= min= max= /autoname;        
    run;
        
    * summarize ctpp auto; 
    proc means data=OUTDATA.proxmeas_national 
        (where=(o_stfid ne d_stfid and flow_ctpp_auto=1)) n median mean std min max;
        title "summarize ctpp auto";
        weight jobs;
        var dist auto speed_auto ctpp_auto speed_ctpp_auto
            ;
        output out=access_means_autoctpp n= median= mean= std= min= max= /autoname;        
    run;
        
    * summarize ctpp transit; 
    proc means data=OUTDATA.proxmeas_national 
        (where=(o_stfid ne d_stfid and flow_ctpp_tran=1)) n median mean std min max;
        title "summarize ctpp tran";
        weight jobs;
        var dist tran speed_tran ctpp_tran speed_ctpp_tran
            ;
        output out=access_means_tranctpp n= median= mean= std= min= max= /autoname;        
    run;

    proc print data=access_means_proxall;
    run;

    * table;
    data access_means_proxtable;
        length tabset $11;
        set access_means_autompo (in=autompo 
                keep=dist_median auto_median speed_auto_median
                     dist_mean auto_mean speed_auto_mean
                     dist_stddev auto_stddev speed_auto_stddev
                rename=(auto_median=time_median speed_auto_median=speed_time_median
                        auto_mean=time_mean speed_auto_mean=speed_time_mean
                        auto_stddev=time_stddev speed_auto_stddev=speed_time_stddev
                        )
                )
            access_means_proxall (in=autoall
                keep=dist_median auto_median speed_auto_median
                     dist_mean auto_mean speed_auto_mean
                     dist_stddev auto_stddev speed_auto_stddev
                rename=(auto_median=time_median speed_auto_median=speed_time_median
                        auto_mean=time_mean speed_auto_mean=speed_time_mean
                        auto_stddev=time_stddev speed_auto_stddev=speed_time_stddev
                        )
                )
            access_means_tranmpo (in=tranmpo
                keep=dist_median tran_median speed_tran_median
                     dist_mean tran_mean speed_tran_mean
                     dist_stddev tran_stddev speed_tran_stddev
                rename=(tran_median=time_median speed_tran_median=speed_time_median
                        tran_mean=time_mean speed_tran_mean=speed_time_mean
                        tran_stddev=time_stddev speed_tran_stddev=speed_time_stddev
                        )
                )
            access_means_proxall (in=tranall
                keep=dist_median tran_median speed_tran_median
                     dist_mean tran_mean speed_tran_mean
                     dist_stddev tran_stddev speed_tran_stddev
                rename=(tran_median=time_median speed_tran_median=speed_time_median
                        tran_mean=time_mean speed_tran_mean=speed_time_mean
                        tran_stddev=time_stddev speed_tran_stddev=speed_time_stddev
                        )
                )
            access_means_autoctpp (in=autoctppmpo
                keep=dist_median auto_median speed_auto_median
                     dist_mean auto_mean speed_auto_mean
                     dist_stddev auto_stddev speed_auto_stddev
                rename=(auto_median=time_median speed_auto_median=speed_time_median
                        auto_mean=time_mean speed_auto_mean=speed_time_mean
                        auto_stddev=time_stddev speed_auto_stddev=speed_time_stddev
                        )
                )
            access_means_autoctpp (in=autoctppsur
                keep=Dist_median ctpp_auto_median speed_ctpp_auto_median
                     dist_mean ctpp_auto_mean speed_ctpp_auto_mean
                     dist_stddev ctpp_auto_stddev speed_ctpp_auto_stddev
                rename=(ctpp_auto_median=time_median speed_ctpp_auto_median=speed_time_median
                        ctpp_auto_mean=time_mean speed_ctpp_auto_mean=speed_time_mean
                        ctpp_auto_stddev=time_stddev speed_ctpp_auto_stddev=speed_time_stddev
                        )
                )
            access_means_tranctpp (in=tranctppmpo
                keep=dist_median tran_median speed_tran_median
                     dist_mean tran_mean speed_tran_mean
                     dist_stddev tran_stddev speed_tran_stddev
                rename=(tran_median=time_median speed_tran_median=speed_time_median
                        tran_mean=time_mean speed_tran_mean=speed_time_mean
                        tran_stddev=time_stddev speed_tran_stddev=speed_time_stddev
                        )
                )
            access_means_tranctpp (in=tranctppsur
                keep=Dist_median ctpp_tran_median speed_ctpp_tran_median
                     dist_mean ctpp_tran_mean speed_ctpp_tran_mean
                     dist_stddev ctpp_tran_stddev speed_ctpp_tran_stddev
                rename=(ctpp_tran_median=time_median speed_ctpp_tran_median=speed_time_median
                        ctpp_tran_mean=time_mean speed_ctpp_tran_mean=speed_time_mean
                        ctpp_tran_stddev=time_stddev speed_ctpp_tran_stddev=speed_time_stddev
                        )
                )
            ;
        tabset='none';
        if (autompo) then tabset='autompo'; * Table C1, Column 1;
        if (autoall) then tabset='autoall'; * Table C1, Column 2;
        if (tranmpo) then tabset='tranmpo'; * Table C1, Column 5;
        if (tranall) then tabset='tranall'; * Table C1, Column 6;
        if (autoctppmpo) then tabset='autoctppmpo'; * Table C1, Column 3;
        if (autoctppsur) then tabset='autoctppsur'; * Table C1, Column 4;
        if (tranctppmpo) then tabset='tranctppmpo'; * Table C1, Column 7;
        if (tranctppsur) then tabset='tranctppsur'; * Table C1, Column 8;
    run;
    proc print data=access_means_proxtable;
        title 'Table C1: summary of travel times';
    run;
            
    * summarize auto non-impute; 
    proc means data=OUTDATA.proxmeas_national 
        (where=(o_stfid ne d_stfid and predict_auto=0)) n sum;
        title "summarize auto non-impute";
        var jobs
            ;
        output out=access_means n= sum= /autoname;        
    run;
        
    * summarize transit non-impute; 
    proc means data=OUTDATA.proxmeas_national 
        (where=(o_stfid ne d_stfid and predict_tran=0 and speed_tran gt 3)) n sum;
        title "summarize transit non-impute";
        var jobs
            ;
        output out=access_means n= sum= /autoname;        
    run;
                        
    * summarize all tracts with flows; 
    proc means data=OUTDATA.proxmeas_national 
        (where=(o_stfid ne d_stfid)) n sum;
        title "summarize all tracts with flows";
        var jobs
            ;
        output out=access_means n= sum= /autoname;        
    run;
        
    * summarize ctpp auto; 
    proc means data=OUTDATA.proxmeas_national 
        (where=(o_stfid ne d_stfid and flow_ctpp_auto=1)) n sum;
        title "summarize ctpp auto";
        var jobs
            ;
        output out=access_means n= sum= /autoname;        
    run;
        
    * summarize ctpp transit; 
    proc means data=OUTDATA.proxmeas_national 
        (where=(o_stfid ne d_stfid and flow_ctpp_tran=1)) n sum;
        title "summarize ctpp tran";
        var jobs
            ;
        output out=access_means n= sum= /autoname;        
    run;
        
    * summarize all tracts with flows out of tract; 
    proc means data=OUTDATA.proxmeas_national 
        (where=(o_stfid ne d_stfid)) n sum;
        title "summarize all tracts with flows out of tract";
        var jobs
            ;
        output out=access_means n= sum= /autoname;        
    run;
        
    * summarize all tracts with flows in tract; 
    proc means data=OUTDATA.proxmeas_national 
        (where=(o_stfid=d_stfid)) n sum;
        title "summarize all tracts with flows in tract";
        var jobs
            ;
        output out=access_means n= sum= /autoname;        
    run;
        
    * estimation of ctpp travel time based on mpo times - auto;
    proc reg data=OUTDATA.proxmeas_national 
        (where=(o_stfid ne d_stfid and flow_ctpp_auto=1 and predict_auto=0));
        title "regression of auto, ctpp travel time for national";
        weight jobs;
        model ctpp_auto=auto;
    run;
    proc corr data=OUTDATA.proxmeas_national
        (where=(o_stfid ne d_stfid and flow_ctpp_auto=1 and predict_auto=0));
        title "correlation of auto, ctpp travel time for national";
        var ctpp_auto auto;
    run;
    * estimation of ctpp travel time based on mpo times - transit;
    proc reg data=OUTDATA.proxmeas_national 
        (where=(o_stfid ne d_stfid and flow_ctpp_tran=1 and predict_tran=0));
        title "regression of tran, ctpp travel time for national";
        weight jobs;
        model ctpp_tran=tran;
    run;
    proc corr data=OUTDATA.proxmeas_national
        (where=(o_stfid ne d_stfid and flow_ctpp_tran=1 and predict_tran=0));
        title "correlation of tran, ctpp travel time for national";
        var ctpp_tran tran;
    run;

%end;

%mend us_access;

